# 05_Urban_Rural.R
# Created: 2021-12-14 Taka-aki Asano
# Last Modified: 2023-07-29

# package
require("readr")
require("dplyr")
require("tidyr")
require("ggplot2")
require("ggpubr")
require("cjoint")
require("cregg")
theme_set(theme_minimal(base_size = 12))


# load dataset
data <- read_csv("HLWconjoint20210330.csv")


# processing of data
data2 <- select(data, 
                no, time, rl, group, age, psm01, psm02, psm03, psm04, 
                psm05, psm06, psm07, psm08, psm09, psm10, psm11, 
                psm12, psm13, psm14, psm15, psm16, hlw_su, hlw_fsu, 
                hlw_dis, hlw_pol, hlw_so, hlw_ad, hlw_pop, hlw_inc)
data2$hlw_fsu <- ifelse(data2$hlw_fsu == data2$rl, 1, 0)
data2 <- mutate(data2, 
                `Distance` = recode(hlw_dis, `1` = "less than 1 km", 
                                    `2` = "1-5 km", `3` = "5-20 km", 
                                    `4` = "20 km or more"), 
                `Politics` = recode(hlw_pol, `1` = "job creation", 
                                    `2` = "public services", 
                                    `3` = "local economy", 
                                    `4` = "society"), 
                `Social` = recode(hlw_so, `1` = "70% agreed, 30% disagreed", 
                                  `2` = "50% agreed, 50% disagreed", 
                                  `3` = "30% agreed, 70% disagreed"), 
                `Administration` = recode(hlw_ad, `1` = "none", 
                                          `2` = "childbirth and childcare", 
                                          `3` = "elderly care", `4` = "education", 
                                          `5` = "public facility", `6` = "water charge"), 
                `Population` = recode(hlw_pop, `1` = "19,000", 
                                      `2` = "20,000", `3` = "21,000"), 
                `Income` = recode(hlw_inc, `1` = "3.8 million yen", 
                                  `2` = "4 million yen", `3` = "4.2 million yen"))
data2$Distance <- factor(
  data2$Distance, 
  levels = c("less than 1 km", "1-5 km", "5-20 km", "20 km or more")
)
data2$Politics <- factor(
  data2$Politics, 
  levels = c("society", "job creation", "public services", "local economy")
)
data2$Social <- factor(
  data2$Social, 
  levels = c("30% agreed, 70% disagreed", "50% agreed, 50% disagreed", "70% agreed, 30% disagreed")
)
data2$Administration <- factor(
  data2$Administration, 
  levels = c("none", "childbirth and childcare", "elderly care", 
             "education", "public facility", "water charge")
)
data2$Population <- factor(
  data2$Population, 
  levels = c("19,000", "20,000", "21,000")
)
data2$Income <- factor(
  data2$Income, 
  levels = c("3.8 million yen", "4 million yen", "4.2 million yen")
)


# Urban - Rural
data2$Residence <- ifelse(data2$group %in% 1:3, "Urban", "Rural")
data2$Residence <- factor(data2$Residence, levels = c("Rural", "Urban"))

## forced-choice, AMCEs
res1 <- cj(
  data2, 
  hlw_fsu ~ Distance + Politics + Social + Administration + Population + Income, 
  id = ~ no, estimate = "amce", by = ~ Residence, level_order = "descending"
)
res1$level <- as.character(res1$level)
res1 <- rbind(
  res1, 
  c("Rural", "hlw_fsu", "amce", "Distance", "(Distance attribute)", NA, NA, NA, NA, NA, NA, "Rural"), 
  c("Rural", "hlw_fsu", "amce", "Politics", "(Political attribute)", NA, NA, NA, NA, NA, NA, "Rural"), 
  c("Rural", "hlw_fsu", "amce", "Social", "(Social attribute)", NA, NA, NA, NA, NA, NA, "Rural"), 
  c("Rural", "hlw_fsu", "amce", "Administration", "(Administrative attribute)", NA, NA, NA, NA, NA, NA, "Rural"), 
  c("Rural", "hlw_fsu", "amce", "Population", "(Population attribute)", NA, NA, NA, NA, NA, NA, "Rural"), 
  c("Rural", "hlw_fsu", "amce", "Income", "(Average income attribute)", NA, NA, NA, NA, NA, NA, "Rural")
)
res1$level <- factor(
  res1$level, 
  levels = c(
    "20 km or more", "5-20 km", "1-5 km", "less than 1 km", "(Distance attribute)", 
    "local economy", "public services", "job creation", "society", "(Political attribute)", 
    "70% agreed, 30% disagreed", "50% agreed, 50% disagreed", "30% agreed, 70% disagreed", "(Social attribute)", 
    "water charge", "public facility", "education", "elderly care", "childbirth and childcare", "none", "(Administrative attribute)",  
    "21,000", "20,000", "19,000", "(Population attribute)", 
    "4.2 million yen", "4 million yen", "3.8 million yen", "(Average income attribute)"
  )
)
res1$estimate <- as.numeric(res1$estimate)
res1$std.error <- as.numeric(res1$std.error)
res1$z <- as.numeric(res1$z)
res1$p <- as.numeric(res1$p)
res1$lower <- as.numeric(res1$lower)
res1$upper <- as.numeric(res1$upper)
res1.plot <- ggplot(
  res1, aes(x = reorder(level, as.numeric(feature) * -1), y = estimate, 
            ymax = upper, ymin = lower, shape = feature, color = Residence)) + 
  geom_hline(yintercept = 0, color = "gray60") + 
  geom_pointrange(position = position_dodge(width = 0.7), size = 0.2) + 
  scale_colour_manual(values = c("black", "grey60"), na.translate = FALSE) + 
  guides(shape = "none") + coord_flip() + ylim(-0.4, 0.4) + 
  labs(y = "Estimated AMCE") + 
  theme_minimal(base_size = 12) + 
  theme(plot.title = element_blank(), 
        axis.title.y = element_blank(), 
        legend.position = "bottom", 
        legend.title = element_blank())
plot(res1.plot)

diff1 <- cj(
  data2, 
  hlw_fsu ~ Distance + Politics + Social + Administration + Population + Income, 
  id = ~ no, estimate = "amce_diff", by = ~ Residence, level_order = "descending"
)
diff1$level <- as.character(diff1$level)
diff1 <- rbind(
  diff1, 
  c("Urban - Rural", "hlw_fsu", "amce", "Distance", "(Distance attribute)", "Urban", NA, NA, NA, NA, NA, NA), 
  c("Urban - Rural", "hlw_fsu", "amce", "Politics", "(Political attribute)", "Urban", NA, NA, NA, NA, NA, NA), 
  c("Urban - Rural", "hlw_fsu", "amce", "Social", "(Social attribute)", "Urban", NA, NA, NA, NA, NA, NA), 
  c("Urban - Rural", "hlw_fsu", "amce", "Administration", "(Administrative attribute)", "Urban", NA, NA, NA, NA, NA, NA), 
  c("Urban - Rural", "hlw_fsu", "amce", "Population", "(Population attribute)", "Urban", NA, NA, NA, NA, NA, NA), 
  c("Urban - Rural", "hlw_fsu", "amce", "Income", "(Average income attribute)", "Urban", NA, NA, NA, NA, NA, NA)
)
diff1$level <- factor(
  diff1$level, 
  levels = c(
    "20 km or more", "5-20 km", "1-5 km", "less than 1 km", "(Distance attribute)", 
    "local economy", "public services", "job creation", "society", "(Political attribute)", 
    "70% agreed, 30% disagreed", "50% agreed, 50% disagreed", "30% agreed, 70% disagreed", "(Social attribute)", 
    "water charge", "public facility", "education", "elderly care", "childbirth and childcare", "none", "(Administrative attribute)",  
    "21,000", "20,000", "19,000", "(Population attribute)", 
    "4.2 million yen", "4 million yen", "3.8 million yen", "(Average income attribute)"
  )
)
diff1$estimate <- as.numeric(diff1$estimate)
diff1$std.error <- as.numeric(diff1$std.error)
diff1$z <- as.numeric(diff1$z)
diff1$p <- as.numeric(diff1$p)
diff1$lower <- as.numeric(diff1$lower)
diff1$upper <- as.numeric(diff1$upper)
diff1.plot <- ggplot(
  diff1, aes(x = reorder(level, as.numeric(feature) * -1), y = estimate, 
             ymax = upper, ymin = lower, shape = feature)) + 
  geom_hline(yintercept = 0, color = "gray60") + 
  geom_pointrange(size = 0.2) + coord_flip() + 
  ylim(-0.25, 0.25) + facet_wrap(. ~ BY, ncol = 1) + 
  labs(y = "Estimated Difference in AMCEs") + 
  theme_minimal(base_size = 12) + 
  theme(axis.title.y = element_blank(), 
        legend.position = "none", 
        legend.title = element_blank())
plot(diff1.plot)

## forced-choice, MMs
res2 <- cj(
  data2, 
  hlw_fsu ~ Distance + Politics + Social + Administration + Population + Income, 
  id = ~ no, estimate = "mm", by = ~ Residence, level_order = "descending"
)
res2$level <- as.character(res2$level)
res2 <- rbind(
  res2, 
  c("Urban", "hlw_fsu", "mm", "Distance", "(Distance attribute)", NA, NA, NA, NA, NA, NA, "Urban"), 
  c("Urban", "hlw_fsu", "mm", "Politics", "(Political attribute)", NA, NA, NA, NA, NA, NA, "Urban"), 
  c("Urban", "hlw_fsu", "mm", "Social", "(Social attribute)", NA, NA, NA, NA, NA, NA, "Urban"), 
  c("Urban", "hlw_fsu", "mm", "Administration", "(Administrative attribute)", NA, NA, NA, NA, NA, NA, "Urban"), 
  c("Urban", "hlw_fsu", "mm", "Population", "(Population attribute)", NA, NA, NA, NA, NA, NA, "Urban"), 
  c("Urban", "hlw_fsu", "mm", "Income", "(Average income attribute)", NA, NA, NA, NA, NA, NA, "Urban")
)
res2$level <- factor(
  res2$level, 
  levels = c(
    "20 km or more", "5-20 km", "1-5 km", "less than 1 km", "(Distance attribute)", 
    "local economy", "public services", "job creation", "society", "(Political attribute)", 
    "70% agreed, 30% disagreed", "50% agreed, 50% disagreed", "30% agreed, 70% disagreed", "(Social attribute)", 
    "water charge", "public facility", "education", "elderly care", "childbirth and childcare", "none", "(Administrative attribute)",  
    "21,000", "20,000", "19,000", "(Population attribute)", 
    "4.2 million yen", "4 million yen", "3.8 million yen", "(Average income attribute)"
  )
)
res2$estimate <- as.numeric(res2$estimate)
res2$std.error <- as.numeric(res2$std.error)
res2$z <- as.numeric(res2$z)
res2$p <- as.numeric(res2$p)
res2$lower <- as.numeric(res2$lower)
res2$upper <- as.numeric(res2$upper)
res2.plot <- ggplot(
  res2, aes(x = reorder(level, as.numeric(feature) * -1), y = estimate, 
            ymax = upper, ymin = lower, shape = feature, color = Residence)) + 
  geom_hline(yintercept = 0.5, color = "gray60") + 
  geom_pointrange(position = position_dodge(width = 0.7), size = 0.2) + 
  scale_colour_manual(values = c("black", "grey60"), na.translate = FALSE) + 
  guides(shape = "none") + coord_flip() + ylim(0.25, 0.75) + 
  labs(y = "Estimated Marginal Mean") + 
  theme_minimal(base_size = 12) + 
  theme(plot.title = element_blank(), 
        axis.title.y = element_blank(), 
        legend.position = "bottom", 
        legend.title = element_blank())
plot(res2.plot)

diff2 <- cj(
  data2, 
  hlw_fsu ~ Distance + Politics + Social + Administration + Population + Income, 
  id = ~ no, estimate = "mm_diff", by = ~ Residence, level_order = "descending"
)
diff2$level <- as.character(diff2$level)
diff2 <- rbind(
  diff2, 
  c("Urban - Rural", "hlw_fsu", "mm", "Distance", "(Distance attribute)", NA, NA, NA, NA, NA, NA, "Urban"), 
  c("Urban - Rural", "hlw_fsu", "mm", "Politics", "(Political attribute)", NA, NA, NA, NA, NA, NA, "Urban"), 
  c("Urban - Rural", "hlw_fsu", "mm", "Social", "(Social attribute)", NA, NA, NA, NA, NA, NA, "Urban"), 
  c("Urban - Rural", "hlw_fsu", "mm", "Administration", "(Administrative attribute)", NA, NA, NA, NA, NA, NA, "Urban"), 
  c("Urban - Rural", "hlw_fsu", "mm", "Population", "(Population attribute)", NA, NA, NA, NA, NA, NA, "Urban"), 
  c("Urban - Rural", "hlw_fsu", "mm", "Income", "(Average income attribute)", NA, NA, NA, NA, NA, NA, "Urban")
)
diff2$level <- factor(
  diff2$level, 
  levels = c(
    "20 km or more", "5-20 km", "1-5 km", "less than 1 km", "(Distance attribute)", 
    "local economy", "public services", "job creation", "society", "(Political attribute)", 
    "70% agreed, 30% disagreed", "50% agreed, 50% disagreed", "30% agreed, 70% disagreed", "(Social attribute)", 
    "water charge", "public facility", "education", "elderly care", "childbirth and childcare", "none", "(Administrative attribute)",  
    "21,000", "20,000", "19,000", "(Population attribute)", 
    "4.2 million yen", "4 million yen", "3.8 million yen", "(Average income attribute)"
  )
)
diff2$estimate <- as.numeric(diff2$estimate)
diff2$std.error <- as.numeric(diff2$std.error)
diff2$z <- as.numeric(diff2$z)
diff2$p <- as.numeric(diff2$p)
diff2$lower <- as.numeric(diff2$lower)
diff2$upper <- as.numeric(diff2$upper)
diff2.plot <- ggplot(
  diff2, aes(x = reorder(level, as.numeric(feature) * -1), y = estimate, 
             ymax = upper, ymin = lower, shape = feature)) + 
  geom_hline(yintercept = 0, color = "gray60") + 
  geom_pointrange(size = 0.2) + coord_flip() + 
  ylim(-0.15, 0.15) + facet_wrap(. ~ BY, ncol = 1) + 
  labs(y = "Estimated Difference in Marginal Means") + 
  theme_minimal(base_size = 12) + 
  theme(axis.title.y = element_blank(), 
        legend.position = "none", 
        legend.title = element_blank())
plot(diff2.plot)
